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We introduce a new class of quantum Monte Carlo methods, based on a Gaussian quantum opera- 
tor representation of fermionic states. The methods enable first-principles dynamical or equilibrium 
calculations in many-body Fermi systems, and, combined with the existing Gaussian representa- 
tion for bosons, provide a unified method of simulating Bose-Fermi systems. As an application, we 
calculate finite-temperature properties of the two dimensional Hubbard model. 
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Calculating the quantum many-body physics of inter- 
acting Fermi systems is one of the great challenges in 
modern theoretical physics. These issues appear in phys- 
ical problems at all energy scales, from ultra-cold atomic 
physics to high-energy lattice QCD. In even the simplest 
cases, first-principles calculations are inhibited by the 
complexity of the fermionic wavefunction, manifest no- 
toriously in the Fermi sign problem. In previous quan- 
tum Monte Carlo (QMC) techniques, the sign problem 
appears as trajectories with negative weights, which con- 
tribute to a large sampling error pj]. QMC methods are 
also complicated by the calculation of large determinants. 

In this letter, we introduce a new QMC method 
for simulating many-body fermion systems, based on a 
Gaussian phase-space representation. As an application 
to condensed matter and AMO physics, we study the 
well-known Hubbard model. Although it is the sim- 
plest model of interacting fermions on a lattice, it is 
rich in physics and may even describe high-temperature 
superconductivity 0- We show that for the Hubbard 
model the Gaussian representation leads to imaginary- 
time equations with no negative probabilities or weights. 
We demonstrate that this removes the well-known Fermi 
sign problem P, 0,0, by first principles numerical simula- 
tion without fixed-node |j| or variational approximations. 

Phase-space methods |5J provide a way to simulate 
quantum many-body systems both dynamically and at 
finite temperature, and have proved useful in bosonic 
cases. These methods sample the time evolution of a 
positive distribution on an overcomplete basis set, which 
is usually the set of coherent states. However, whereas 
coherent state representations are well-defined in the 
bosonic case, the only known coherent state techniques 
for fermions involve Grassmann algebra^ , which has an 
enormous computational complexity. 

Here we introduce a phase-space method that over- 
comes the problem of Grassmann complexity, using a 
Gaussian expansion for fermions. The operator basis is 
constructed from pairs of Fermi operators. Because these 
pairs obey commutation relations, a natural solution of 
the Grassmann problem is achieved. Furthermore, the 
resulting equations obviate the need to evaluate large de- 
terminants. The elimination of anti-commutators means 
that the technique is far more efficient than previous 
QMC and stochastic fermion methods 7] . We give ex- 
amples in cases of experimental relevance involving the 



dynamical problem of Pauli blocking in molecular disso- 
ciation, and finite temperature correlations of fermions 
in an optical lattice, where the results agree with those 
of other exact methods. We also perform larger simu- 
lations of the 2D Hubbard model, in cases where severe 
sign problems were found previously. 

Our starting point is a general expansion of the system 
density operator: 



P{X,t)A(X)d\ 



(1) 



where P( A , t) is a probability distribution, A is a suitable 
basis for the class of density matrices being considered, 
and d A is the integration measure for the corresponding 
generalized phase-space coordinate A . The operators A 
are non-Hermitian and form a complete basis for density 
operator. Existing phase-space methods are defined for 
systems of bosons, with A constructed of bosonic ladder 
operators^. 

To achieve a unified representation, we define the op- 
erator basis A to be the product of Gaussian forms of 
bosonic and fermionic creation and annihilation opera- 
tors: A = f2A(,A/, where A& and Af are Gaussian forms 
over Mb bosonic modes and M fermionic modes, respec- 
tively, and where the (possibly) complex number Q is 
an additional weighting factor. The properties of the 
bosonic Gaussian representation are given in ||. Here 
we summarise the relevant properties of the Fermionic 
Gaussian form; explicit proofs will be given elsewhere. 

For a system that can be decomposed into M single- 
particle modes, we define a as a column vector of the M 
annihilation operators, and at as the corresponding row 
vector of M creation operators, whose anticommutation 



relations: [afc, aj] + = Skj ■ We also introduce an extended 

(a, (at) T ) , with 



a 



2M-vector of all the operators: 

adjoint defined as a? = (at , a T ) . A general, normally 
ordered Gaussian operator can then be written 



A/ = Pf 



OA 



exp 



(2) 



which, because it is constructed from pairs of operators, 
contain no Grassmann variables. The normalisation, cho- 
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sen to ensure that Tr A/ = 1 , consists of the Pfafhan of 
an antisymmetric form a a of the covariance^. 

Normal ordering, denoted by is defined as in 

the bosonic case, with all annihilation operators to the 
right of the creation operators, except that each pairwise 
reordering involved induces a sign change, e. g. : aia^ : 
= — a/jCLi . We define aniinormal ordering similarly, and 

denote it via curly braces: {aja^} = — a^a] . More gen- 
erally, we can define nested orderings, in which the outer 
ordering does not reorder the inner one. For example, 
{: Acij :a,i} = — OiSjA , where we assume that the kernel 

A always remains normally ordered. 

The generalized covariance a and constant matrix I 
are 2M x 2M matrices, which we can write as 



I n 2 
m 



rn 
n I 



I = 



I 

—I 



(3) 



where the number correlation n is a complex M x M 
matrix, the squeezing correlations m, m + are two inde- 
pendent antisymmetric complex M x M matrices, and I 
is the M-mode identity matrix. 

The phase space of the fermionic representation is A = 
(fi, n, m, m + ), which has a dimension of 1 + p = 1 + 
M(2M — 1) . For a combined Bose-Fermi system, there 
will be an additional M{,(2M{, + 3) bosonic dimensions. 

Under the Gaussian representation, physical quanti- 
ties (operator expectation values) appear as moments of 
the (weighted) distribution ttP, denoted as (..) P . For 
quadratic products, 
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(4) 



For higher order products, the corresponding moments 
can be determined by evaluation of the appropriate 
(Grassmann) Gaussian integral. Note that there is no 
way to calculate the expectation value of single ladder 
operators (or any product that is of odd order) ; the Gaus- 
sian form cannot represent density operators that con- 
tain an odd number of operators. However, in physical 
Hamiltonians, Fermi operators appear only in pairs, and 
so such 'odd' states will not be generated in the course 
of the evolution. For other, physical states, the Gaussian 
basis provides an (over)complete representation. 

To enable canonical or dynamical simulations, we need 
identities that describe the action of operators on the 
density operator as derivatives on elements of the Gaus- 
sian basis. With our ordering notation given above, all 
of the necessary operator identities can be written in a 
compact matrix form: 



A 



{a:o f A:} 
{aa f A} 



— a A + a-— a , 
— —ag_— 

aA-(a- I) — a 



dA 



(g-l)A+(g-l) — (^-I).(5) 



The matrix derivative is here defined as (d/da)^ — 
djduu^ . Notice that there are no determinants (or Pfaf- 
fians) to be calculated in these identities. 

As an application of the fermionic representation, con- 
sider the well-known Hubbard model: H (n-| , nj ) = 



v t,J.<T 



(6) 



where rii,j,v = a| )(T Oj- lCr ={n <r }y , t is the hopping, or 
tunclling, strength, U is strength of on-site interacta- 
tions and /i is the chemical potential, included to con- 
trol the total particle number. The index a denotes spin 
(1 , J ) and the indices i,j label lattice location, with 
denoting a sum over nearest neighbours. The Hubbard 
model is the simplest nontrivial model for strongly in- 
teracting electrons and is thus an important system in 
condensed matter physics, with relevance to the theory 
of high-temperature superconductors 0- It also describes 
an ultracold fermi gas in a optical lattice potential. The 
physics of the model is not yet fully understood , an d al- 
though there are known solutions in the ID case[l(J, this 
is not so for higher dimensions. 

The 2D problem in particular is an important test- 
ing ground for QMC methods. Traditional methods are 
prone to sign problems in the repulsive case (U > 0) 
away from half filling. These are particularly severe for 
large systems, higher dimensio ns, stronger interaction 
and open-shell configurations pi, ITlj . 

The equilibrium state at temperature T — 1/kp.T can 
be cast into an 'imaginary time' differential equation for 
the unnormalised density operator: dp/dr — — ^[H ,p\+ ■ 
We make use of the representation by expanding the den- 
sity operator in terms of the Gaussian operators and 
applying the identities in Eq. (jjjj. After integrating 
by parts, we arrive at an equation for the distribution 
function, which we can sample numerically, by solv- 
ing stochastic phase-space equations. Although there 
is never any need to calculate determinants with these 
methods, the sampling error typically grows in (imagi- 
nary) time unless a suitable choice of 'stochastic gauge' is 
made[l2j. in which one exploits the overcomplete nature 
of the basis to keep the distribution compact. Stochas- 
tic gauges can also be used to eliminate boundary terms 
that may arise in the partial integration step. 

Before applying this procedure to the Hubbard model, 
we first rewrite the interaction terms as 



U : n 



3,3,1 n 3,3,\ 



-{U^-.^-sn^f :, (7) 
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where s = U/\U\. The extra terms here vanish because 
of the anticommuting property of fermion operators, but 
they do lead to additional stochastic terms. Such terms 
are examples of a new type of stochastic gauge, and one 
that is unique to fermions: vanishing operator products 
can be used to modify the stochastic behaviour of the 
phase-space equations without affecting the averaged re- 
sults. With this choice of terms, we map the imaginary- 
time calculation onto a set of real Stratonovich stochastic 
equations, which, in matrix form, are 

^ = i{(I-n £r )A( 1 W + n CT A( 2 )(I-n (T )}.(8) 

Here we have introduced the matrix: A; I „ = 

tf>{i,i) ~ l\U\(sn jd -„ - n jijta + ~) - M + /£j r) j > 

where Suj\ = 1 if the i,j correspond to nearest neigh- 
bour sites and is otherwise 0, and where / = — s for a =| 
and 1 otherwise. The real Gaussian noise if{r) is de- 
fined by the correlations ^/'V')) = 2\U\S{t - 

T ')Sj,j'Sr,r' ■ The weights for each trajectory evolve as 
physically expected for energy-weighted averages, with 
dn/dr = — Sli7(ni , nj ). Because the equations for the 
phase-space variables rii d ^ are all real, the weights will 
all remain positive, thereby avoiding the traditional man- 
ifestation of the sign problem. 

Consider first the case where t — 0, which describes, for 
example, an ultracold Fermi gas in a deep optical lattice 
potential. In Fig. ^ we plot the the correlation function 
= (: nin 2 :) / (n\) (n 2 ) for U > 0, revealing a strong 
antibunching effect at low temperatures. The sampling 
error here is very small, because of the restricted phase- 
space explored by the simulation. 

Whether the method can overcome the fundamental 
cause of the sign problem, which is the complexity of 
fcrmionic states, must be demonstrated by calculating 
physical quantities in cases where the sign problem is 
known to occur in other methods. Thus we calculate the 
total energy for t = 1, U = 4 in two dimensions as a func- 
tion of temperature, for different fillings. The results for a 
16x16 lattice are shown in Fig. [21 in which to obtain good 
sampling with the spreading weights, we use a branching 
technique^j|. For a 4x4 lattice at an inverse tempera- 
ture of t — 7, we calculate E = —13.1 ± 1.2 at n = 0.5 
and E = -19.62 ± 0.87 at n = .313 ± 0.005, which can 
be compared to zero-temperature, exact-diagonalisation 
results, namely E — —13.62 for n = 0.5 (half filling) 
and E = -19.57 for n = 0.3125 (10 atoms) 14] . At a 
filling of n = 0.412 ± 0.01, for which the sign deterio- 
rates for a projector QMC calculation^]], we calculate 
E = — 16.5 ± 1.5. We emphasise that, unlike projector 
QMC |2j , the Gaussian method can calculate any physi- 
cal correlation function, at any temperature. 

As an application of the method to a dynamical calcu- 
lation, in particular to a composite Bose/Fermi system, 
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Figure 1: Second-order correlation function versus in- 
verse temperature r for t = 0, U = 2 and /i = 1, for which 
(rij) = 0.5. The solid curve gives the simulation result, and 
the dashed and dot-dashed line show the estimated sampling 
error and deviation from the analytic result, respectively (on 
a xlOOO scale). Calculated from 100,000 trajectories 
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Figure 2: Total energy E per site versus inverse temperature r 
for a 16 x 16 2D lattice for chemical potentials /i = 2 (solid), 
fi — 1 (dashed) and fi — (dot-dashed). Curves without 
crosses give the number of particles per site for /i = 1 (dashed) 
and = (dot-dashed). U = 4, t = 1, and 50 paths initially. 
Dotted curves give an estimate of sampling error. 



we consider the process of the dissociation of a molec- 
ular Bose condensate into its constituent atoms, which 
may be fermions or bosons. For simplicity, we consider 
two atomic modes, representing, for example, states of 
different spin or momenta, coupled to a single molecular 
mode, via the effective interaction H = d^bib2 + h.c, 
where foj, bj are the atomic creation and annihilation 
operators and d\ a, are the bosonic molecular opera- 
tors. Realistic models of the atomic-molecular Feshbach 
resonances contain such terms to describe the coupling, 
and it is important to illustrate how this method can 
represent them. Because the normal spin-spin correla- 
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Figure 3: Molecular dissociation into pairs of fermionic (solid 
line) or bosonic (dot-dashed line) atoms. For the fermionic 
case, the dashed curve gives the truncated number-state cal- 
culation, and the dotted lines the estimated sampling error. 
In the bosonic case, the estimated sampling error is too small 
to be distinctly plotted on this graph. The initial state is 
a molecular coherent state ( jV mo iecuies(0) = 9). Calculated 
from 10,000 trajectories. 

tions < b\i)2 > will remain zero in this system (if ini- 
tially zero), the phase space of the system reduces to 
A = (a, a + , ni,ri2, to, m + ) , i. c. four complex atomic 
variables and two complex molecular amplitudes. Ap- 
plying the identities in Eq. J5J and in Q, we derive the 
following phase-space equations for the time evolution, 
(where H ► bosonic, > fermionic): 

rij = i(a + m — am + ) ± y/inj (to£* + m + Q) , 
to = —ia(l ± n\ ± 77.2) + \~i (±to 2 Cj + nir^CI) > 
m+ = ia + (l ± ni ± n 2 ) + Vi {nm2(,* ± rn^Q) , 
a = —im — V^Ci 7 a+ = * m+ + v^Cs , (9) 
where j — 1,2 and where the Cfc(i) are two complex Gaus- 



sian noises, defined by the correlations (Cfc (^)Cfc' (*')) = 
0, (Ck{t)Q,(t')} =S ktk >5(t-t') . Simulations of Eqs. © 
are compared with truncated number-state-based calcu- 
lations in Fig. [3J Although the initial rates of conver- 
sion are the same in each case, a Pauli blocking effect 
soon slows the fermionic conversion, in contrast to an 
enhanced bosonic conversion. We note that in these real- 
time calculations, a growing sampling error appears to 
be a generic property, although a prudent gauge choice 
may control the growth rate for a certain time. 

In summary, we have introduced here an operator rep- 
resentation that is able to represent arbitrary physical 
states of fermions. Together with the corresponding 
bosonic representation, it is the largest class of repre- 
sentations that can be constructed using an operator ba- 
sis that is Gaussian in the ladder operators. We have 
presented identities for first-principles calculations of the 
time evolution of quantum systems, both dynamical (real 
time) and canonical (imaginary time) . A quadratic mas- 
ter equation maps to deterministic equations, whereas 
interacting systems with quartic terms in the Hamilto- 
nian generate stochastic equations, provided a suitable 
stochastic gauge is chosen that eliminates all boundary 
terms. No computationally intensive determinant calcu- 
lations are involved. 

The simple examples given here show how the one, uni- 
fied method can solve both fermionic and bosonic prob- 
lems, making it well suited to simulating Bose-Fermi mix- 
tures, and to studying from first-principles, for exam- 
ple, the BEC/BCS crossover. Importantly, a new type of 
Fermi stochastic freedom can be used to map canonical 
calculations of the Hubbard type onto a real subspace. 
We have thereby been able to numerically simulate the 
Hubbard model without sign error, even without employ- 
ing any of the sophisticated sampling techniques that 
have been developed over time to optimise more conven- 
tional QMC methods. The application of such techniques 
to the Gaussian approach is yet to be explored. 
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